#include <math.h>
#include "utility.h"

#define MATH_PI 3.141592653589793238462643383279

namespace libvec
{
  float radians(float degrees)
  {
    return degrees / 180.0f * MATH_PI;
  }

  double radians(double degrees)
  {
    return degrees / 180.0f * MATH_PI;
  }

  float degrees(float radians)
  {
    return radians / MATH_PI * 180.0f;
  }

  double degrees(double radians)
  {
    return radians / MATH_PI * 180.0f;
  }

  float length(vec2 v)
  {
    return sqrt(v.x*v.x + v.y*v.y);
  }

  float length(vec3 v)
  {
    return sqrt(v.x*v.x + v.y*v.y + v.z*v.z);
  }

  float length(vec4 v)
  {
    return sqrt(v.x*v.x + v.y*v.y + v.z*v.z + v.w*v.w);
  }

  double length(dvec2 v)
  {
    return sqrt(v.x*v.x + v.y*v.y);
  }

  double length(dvec3 v)
  {
    return sqrt(v.x*v.x + v.y*v.y + v.z*v.z);
  }

  double length(dvec4 v)
  {
    return sqrt(v.x*v.x + v.y*v.y + v.z*v.z + v.w*v.w);
  }

  float distance(vec2 p0, vec2 p1)
  {
    return length(p0 - p1);
  }

  float distance(vec3 p0, vec3 p1)
  {
    return length(p0 - p1);
  }

  float distance(vec4 p0, vec4 p1)
  {
    return length(p0 - p1);
  }

  double distance(dvec2 p0, dvec2 p1)
  {
    return length(p0 - p1);
  }

  double distance(dvec3 p0, dvec3 p1)
  {
    return length(p0 - p1);
  }

  double distance(dvec4 p0, dvec4 p1)
  {
    return length(p0 - p1);
  }

  float dot(vec2 v1, vec2 v2)
  {
    return v1.x*v2.x + v1.y*v2.y;
  }

  float dot(vec3 v1, vec3 v2)
  {
    return v1.x*v2.x + v1.y*v2.y + v1.z*v2.z;
  }

  float dot(vec4 v1, vec4 v2)
  {
    return v1.x*v2.x + v1.y*v2.y + v1.z*v2.z + v1.w*v2.w;
  }

  double dot(dvec2 v1, dvec2 v2)
  {
    return v1.x*v2.x + v1.y*v2.y;
  }

  double dot(dvec3 v1, dvec3 v2)
  {
    return v1.x*v2.x + v1.y*v2.y + v1.z*v2.z;
  }

  double dot(dvec4 v1, dvec4 v2)
  {
    return v1.x*v2.x + v1.y*v2.y + v1.z*v2.z + v1.w*v2.w;
  }

  vec3 cross(vec3 v1, vec3 v2)
  {
    static vec3 out;

    out.x = v1.y*v2.z - v1.z*v2.y;
    out.y = v1.z*v2.x - v1.x*v2.z;
    out.z = v1.x*v2.y - v1.y*v2.x;

    return out;
  }

  dvec3 cross(dvec3 v1, dvec3 v2)
  {
    static dvec3 out;

    out.x = v1.y*v2.z - v1.z*v2.y;
    out.y = v1.z*v2.x - v1.x*v2.z;
    out.z = v1.x*v2.y - v1.y*v2.x;

    return out;
  }

  vec2 normalize(vec2 v)
  {
    static vec2 out;

    out = v / length(v);

    return out;
  }

  vec3 normalize(vec3 v)
  {
    static vec3 out;

    out = v / length(v);

    return out;
  }

  vec4 normalize(vec4 v)
  {
    static vec4 out;

    out = v / length(v);

    return out;
  }

  dvec2 normalize(dvec2 v)
  {
    static dvec2 out;

    out = v / length(v);

    return out;
  }

  dvec3 normalize(dvec3 v)
  {
    static dvec3 out;

    out = v / length(v);

    return out;
  }

  dvec4 normalize(dvec4 v)
  {
    static dvec4 out;

    out = v / length(v);

    return out;
  }

  vec2 faceforward(vec2 n, vec2 i, vec2 nref)
  {
    static vec2 out;

    if (dot(nref, i) < 0)
    {
      out = n;
    }
    else
    {
      out = -n;
    }

    return out;
  }

  vec3 faceforward(vec3 n, vec3 i, vec3 nref)
  {
    static vec3 out;

    if (dot(nref, i) < 0)
    {
      out = n;
    }
    else
    {
      out = -n;
    }

    return out;
  }

  vec4 faceforward(vec4 n, vec4 i, vec4 nref)
  {
    static vec4 out;

    if (dot(nref, i) < 0)
    {
      out = n;
    }
    else
    {
      out = -n;
    }

    return out;
  }

  dvec2 faceforward(dvec2 n, dvec2 i, dvec2 nref)
  {
    static dvec2 out;

    if (dot(nref, i) < 0)
    {
      out = n;
    }
    else
    {
      out = -n;
    }

    return out;
  }

  dvec3 faceforward(dvec3 n, dvec3 i, dvec3 nref)
  {
    static dvec3 out;

    if (dot(nref, i) < 0)
    {
      out = n;
    }
    else
    {
      out = -n;
    }

    return out;
  }

  dvec4 faceforward(dvec4 n, dvec4 i, dvec4 nref)
  {
    static dvec4 out;

    if (dot(nref, i) < 0)
    {
      out = n;
    }
    else
    {
      out = -n;
    }

    return out;
  }

  vec2 reflect(vec2 i, vec2 n)
  {
    static vec2 out;

    float p = 2.0f * dot(i, n);

    out.x = i.x - p * n.x;
    out.y = i.y - p * n.y;

    return out;
  }

  vec3 reflect(vec3 i, vec3 n)
  {
    static vec3 out;

    float p = 2.0f * dot(i, n);

    out.x = i.x - p * n.x;
    out.y = i.y - p * n.y;
    out.z = i.z - p * n.z;

    return out;
  }

  vec4 reflect(vec4 i, vec4 n)
  {
    static vec4 out;

    float p = 2.0f * dot(i, n);

    out.x = i.x - p * n.x;
    out.y = i.y - p * n.y;
    out.z = i.z - p * n.z;
    out.w = i.w - p * n.w;

    return out;
  }

  dvec2 reflect(dvec2 i, dvec2 n)
  {
    static dvec2 out;

    double p = 2.0f * dot(i, n);

    out.x = i.x - p * n.x;
    out.y = i.y - p * n.y;

    return out;
  }

  dvec3 reflect(dvec3 i, dvec3 n)
  {
    static dvec3 out;

    double p = 2.0f * dot(i, n);

    out.x = i.x - p * n.x;
    out.y = i.y - p * n.y;
    out.z = i.z - p * n.z;

    return out;
  }

  dvec4 reflect(dvec4 i, dvec4 n)
  {
    static dvec4 out;

    double p = 2.0f * dot(i, n);

    out.x = i.x - p * n.x;
    out.y = i.y - p * n.y;
    out.z = i.z - p * n.z;
    out.w = i.w - p * n.w;

    return out;
  }

  vec2 refract(vec2 i, vec2 n, float eta)
  {
    static vec2 out;

    float k = 1.0f - eta * eta * (1.0f - dot(n, i) * dot(n, i));

    if (k < 0.0f)
    {
      return vec2(0.0f);
    }

    out = i * eta - n * (eta * dot(n, i) + sqrt(k));

    return out;
  }

  vec3 refract(vec3 i, vec3 n, float eta)
  {
    static vec3 out;

    float k = 1.0f - eta * eta * (1.0f - dot(n, i) * dot(n, i));

    if (k < 0.0f)
    {
      return vec3(0.0f);
    }

    out = i * eta - n * (eta * dot(n, i) + sqrt(k));

    return out;
  }

  vec4 refract(vec4 i, vec4 n, float eta)
  {
    static vec4 out;

    float k = 1.0f - eta * eta * (1.0f - dot(n, i) * dot(n, i));

    if (k < 0.0f)
    {
      return vec4(0.0f);
    }

    out = i * eta - n * (eta * dot(n, i) + sqrt(k));

    return out;
  }

  dvec2 refract(dvec2 i, dvec2 n, double eta)
  {
    static dvec2 out;

    double k = 1.0f - eta * eta * (1.0f - dot(n, i) * dot(n, i));

    if (k < 0.0f)
    {
      return dvec2(0.0f);
    }

    out = i * eta - n * (eta * dot(n, i) + sqrt(k));

    return out;
  }

  dvec3 refract(dvec3 i, dvec3 n, double eta)
  {
    static dvec3 out;

    double k = 1.0f - eta * eta * (1.0f - dot(n, i) * dot(n, i));

    if (k < 0.0f)
    {
      return dvec3(0.0f);
    }

    out = i * eta - n * (eta * dot(n, i) + sqrt(k));

    return out;
  }

  dvec4 refract(dvec4 i, dvec4 n, double eta)
  {
    static dvec4 out;

    double k = 1.0f - eta * eta * (1.0f - dot(n, i) * dot(n, i));

    if (k < 0.0f)
    {
      return dvec4(0.0f);
    }

    out = i * eta - n * (eta * dot(n, i) + sqrt(k));

    return out;
  }

  mat2 matrixCompMult(mat2 x, mat2 y)
  {
    static mat2 out;

    out.v1.x = x.v1.x * y.v1.x;
    out.v1.y = x.v1.y * y.v1.y;

    out.v2.x = x.v2.x * y.v2.x;
    out.v2.y = x.v2.y * y.v2.y;

    return out;
  }

  mat3 matrixCompMult(mat3 x, mat3 y)
  {
    static mat3 out;

    out.v1.x = x.v1.x * y.v1.x;
    out.v1.y = x.v1.y * y.v1.y;
    out.v1.z = x.v1.z * y.v1.y;

    out.v2.x = x.v2.x * y.v2.x;
    out.v2.y = x.v2.y * y.v2.y;
    out.v2.z = x.v2.z * y.v2.z;

    out.v3.x = x.v3.x * y.v3.x;
    out.v3.y = x.v3.y * y.v3.y;
    out.v3.z = x.v3.z * y.v3.z;

    return out;
  }

  mat4 matrixCompMult(mat4 x, mat4 y)
  {
    static mat4 out;

    out.v1.x = x.v1.x * y.v1.x;
    out.v1.y = x.v1.y * y.v1.y;
    out.v1.z = x.v1.z * y.v1.y;
    out.v1.w = x.v1.w * y.v1.w;

    out.v2.x = x.v2.x * y.v2.x;
    out.v2.y = x.v2.y * y.v2.y;
    out.v2.z = x.v2.z * y.v2.z;
    out.v2.w = x.v2.w * y.v2.w;

    out.v3.x = x.v3.x * y.v3.x;
    out.v3.y = x.v3.y * y.v3.y;
    out.v3.z = x.v3.z * y.v3.z;
    out.v3.w = x.v3.w * y.v3.w;

    out.v4.x = x.v4.x * y.v4.x;
    out.v4.y = x.v4.y * y.v4.y;
    out.v4.z = x.v4.z * y.v4.z;
    out.v4.w = x.v4.w * y.v4.w;

    return out;
  }

  dmat2 matrixCompMult(dmat2 x, dmat2 y)
  {
    static dmat2 out;

    out.v1.x = x.v1.x * y.v1.x;
    out.v1.y = x.v1.y * y.v1.y;

    out.v2.x = x.v2.x * y.v2.x;
    out.v2.y = x.v2.y * y.v2.y;

    return out;
  }

  dmat3 matrixCompMult(dmat3 x, dmat3 y)
  {
    static dmat3 out;

    out.v1.x = x.v1.x * y.v1.x;
    out.v1.y = x.v1.y * y.v1.y;
    out.v1.z = x.v1.z * y.v1.y;

    out.v2.x = x.v2.x * y.v2.x;
    out.v2.y = x.v2.y * y.v2.y;
    out.v2.z = x.v2.z * y.v2.z;

    out.v3.x = x.v3.x * y.v3.x;
    out.v3.y = x.v3.y * y.v3.y;
    out.v3.z = x.v3.z * y.v3.z;

    return out;
  }

  dmat4 matrixCompMult(dmat4 x, dmat4 y)
  {
    static dmat4 out;

    out.v1.x = x.v1.x * y.v1.x;
    out.v1.y = x.v1.y * y.v1.y;
    out.v1.z = x.v1.z * y.v1.y;
    out.v1.w = x.v1.w * y.v1.w;

    out.v2.x = x.v2.x * y.v2.x;
    out.v2.y = x.v2.y * y.v2.y;
    out.v2.z = x.v2.z * y.v2.z;
    out.v2.w = x.v2.w * y.v2.w;

    out.v3.x = x.v3.x * y.v3.x;
    out.v3.y = x.v3.y * y.v3.y;
    out.v3.z = x.v3.z * y.v3.z;
    out.v3.w = x.v3.w * y.v3.w;

    out.v4.x = x.v4.x * y.v4.x;
    out.v4.y = x.v4.y * y.v4.y;
    out.v4.z = x.v4.z * y.v4.z;
    out.v4.w = x.v4.w * y.v4.w;

    return out;
  }

  mat2 outerProduct(vec2 c, vec2 r)
  {
    static mat2 out;

    out.v1.x = c.x * r.x;
    out.v1.y = c.y * r.x;

    out.v2.x = c.x * r.y;
    out.v2.y = c.y * r.y;

    return out;
  }

  mat3 outerProduct(vec3 c, vec3 r)
  {
    static mat3 out;

    out.v1.x = c.x * r.x;
    out.v1.y = c.y * r.x;
    out.v1.z = c.z * r.x;

    out.v2.x = c.x * r.y;
    out.v2.y = c.y * r.y;
    out.v2.z = c.z * r.y;

    out.v3.x = c.x * r.z;
    out.v3.y = c.y * r.z;
    out.v3.z = c.z * r.z;

    return out;
  }

  mat4 outerProduct(vec4 c, vec4 r)
  {
    static mat4 out;

    out.v1.x = c.x * r.x;
    out.v1.y = c.y * r.x;
    out.v1.z = c.z * r.x;
    out.v1.w = c.w * r.x;

    out.v2.x = c.x * r.y;
    out.v2.y = c.y * r.y;
    out.v2.z = c.z * r.y;
    out.v2.w = c.w * r.y;

    out.v3.x = c.x * r.z;
    out.v3.y = c.y * r.z;
    out.v3.z = c.z * r.z;
    out.v3.w = c.w * r.z;

    out.v4.x = c.x * r.w;
    out.v4.y = c.y * r.w;
    out.v4.z = c.z * r.w;
    out.v4.w = c.w * r.w;

    return out;
  }

  dmat2 outerProduct(dvec2 c, dvec2 r)
  {
    static dmat2 out;

    out.v1.x = c.x * r.x;
    out.v1.y = c.y * r.x;

    out.v2.x = c.x * r.y;
    out.v2.y = c.y * r.y;

    return out;
  }

  dmat3 outerProduct(dvec3 c, dvec3 r)
  {
    static dmat3 out;

    out.v1.x = c.x * r.x;
    out.v1.y = c.y * r.x;
    out.v1.z = c.z * r.x;

    out.v2.x = c.x * r.y;
    out.v2.y = c.y * r.y;
    out.v2.z = c.z * r.y;

    out.v3.x = c.x * r.z;
    out.v3.y = c.y * r.z;
    out.v3.z = c.z * r.z;

    return out;
  }

  dmat4 outerProduct(dvec4 c, dvec4 r)
  {
    static dmat4 out;

    out.v1.x = c.x * r.x;
    out.v1.y = c.y * r.x;
    out.v1.z = c.z * r.x;
    out.v1.w = c.w * r.x;

    out.v2.x = c.x * r.y;
    out.v2.y = c.y * r.y;
    out.v2.z = c.z * r.y;
    out.v2.w = c.w * r.y;

    out.v3.x = c.x * r.z;
    out.v3.y = c.y * r.z;
    out.v3.z = c.z * r.z;
    out.v3.w = c.w * r.z;

    out.v4.x = c.x * r.w;
    out.v4.y = c.y * r.w;
    out.v4.z = c.z * r.w;
    out.v4.w = c.w * r.w;

    return out;
  }

  mat2 transpose(mat2 m)
  {
    static mat2 out;

    out.v1.x = m.v1.x;
    out.v2.x = m.v1.y;

    out.v1.y = m.v2.x;
    out.v2.y = m.v2.y;

    return out;
  }

  mat3 transpose(mat3 m)
  {
    static mat3 out;

    out.v1.x = m.v1.x;
    out.v2.x = m.v1.y;
    out.v3.x = m.v1.z;

    out.v1.y = m.v2.x;
    out.v2.y = m.v2.y;
    out.v3.y = m.v2.z;

    out.v1.z = m.v3.x;
    out.v2.z = m.v3.y;
    out.v3.z = m.v3.z;

    return out;
  }

  mat4 transpose(mat4 m)
  {
    static mat4 out;

    out.v1.x = m.v1.x;
    out.v2.x = m.v1.y;
    out.v3.x = m.v1.z;
    out.v4.x = m.v1.w;

    out.v1.y = m.v2.x;
    out.v2.y = m.v2.y;
    out.v3.y = m.v2.z;
    out.v4.y = m.v2.w;

    out.v1.z = m.v3.x;
    out.v2.z = m.v3.y;
    out.v3.z = m.v3.z;
    out.v4.z = m.v3.w;

    out.v1.w = m.v4.x;
    out.v2.w = m.v4.y;
    out.v3.w = m.v4.z;
    out.v4.w = m.v4.w;

    return out;
  }

  dmat2 transpose(dmat2 m)
  {
    static dmat2 out;

    out.v1.x = m.v1.x;
    out.v2.x = m.v1.y;

    out.v1.y = m.v2.x;
    out.v2.y = m.v2.y;

    return out;
  }

  dmat3 transpose(dmat3 m)
  {
    static dmat3 out;

    out.v1.x = m.v1.x;
    out.v2.x = m.v1.y;
    out.v3.x = m.v1.z;

    out.v1.y = m.v2.x;
    out.v2.y = m.v2.y;
    out.v3.y = m.v2.z;

    out.v1.z = m.v3.x;
    out.v2.z = m.v3.y;
    out.v3.z = m.v3.z;

    return out;
  }

  dmat4 transpose(dmat4 m)
  {
    static dmat4 out;

    out.v1.x = m.v1.x;
    out.v2.x = m.v1.y;
    out.v3.x = m.v1.z;
    out.v4.x = m.v1.w;

    out.v1.y = m.v2.x;
    out.v2.y = m.v2.y;
    out.v3.y = m.v2.z;
    out.v4.y = m.v2.w;

    out.v1.z = m.v3.x;
    out.v2.z = m.v3.y;
    out.v3.z = m.v3.z;
    out.v4.z = m.v3.w;

    out.v1.w = m.v4.x;
    out.v2.w = m.v4.y;
    out.v3.w = m.v4.z;
    out.v4.w = m.v4.w;

    return out;
  }

  float determinant(mat2 m)
  {
    return m.v1.x*m.v2.y - m.v2.x*m.v1.y;
  }

  float determinant(mat3 m)
  {
    return
      m.v1.x*(m.v2.y*m.v3.z - m.v3.y*m.v2.z) -
      m.v2.x*(m.v1.y*m.v3.z - m.v3.y*m.v1.z) -
      m.v3.x*(m.v1.y*m.v2.z - m.v2.y*m.v1.z);
  }

  float determinant(mat4 m)
  {
    return
      m.v1.x*(
          m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
          m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
          m.v4.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z)
          ) -
      m.v2.x*(
          m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
          m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
          m.v4.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z)
          ) +
      m.v3.x*(
          m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
          m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
          m.v4.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
          ) -
      m.v4.x*(
          m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) -
          m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) +
          m.v3.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
          );
  }

  double determinant(dmat2 m)
  {
    return m.v1.x*m.v2.y - m.v2.x*m.v1.y;
  }

  double determinant(dmat3 m)
  {
    return
      m.v1.x*(m.v2.y*m.v3.z - m.v3.y*m.v2.z) -
      m.v2.x*(m.v1.y*m.v3.z - m.v3.y*m.v1.z) -
      m.v3.x*(m.v1.y*m.v2.z - m.v2.y*m.v1.z);
  }

  double determinant(dmat4 m)
  {
    return
      m.v1.x*(
          m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
          m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
          m.v4.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z)
          ) -
      m.v2.x*(
          m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
          m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
          m.v4.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z)
          ) +
      m.v3.x*(
          m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
          m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
          m.v4.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
          ) -
      m.v4.x*(
          m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) -
          m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) +
          m.v3.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
          );
  }

  mat2 inverse(mat2 m)
  {
    static mat2 out;

    out.v1.x =  m.v2.y/(m.v1.x*m.v2.y - m.v1.y*m.v2.x);
    out.v1.y = -m.v1.y/(m.v1.x*m.v2.y - m.v1.y*m.v2.x);
    out.v2.x = -m.v2.x/(m.v1.x*m.v2.y - m.v1.y*m.v2.x);
    out.v2.y =  m.v1.x/(m.v1.x*m.v2.y - m.v1.y*m.v2.x);

    return out;
  }

  mat3 inverse(mat3 m)
  {
    static mat3 out;

    out.v1.x = (m.v2.y*m.v3.z - m.v2.z*m.v3.y) / (
        m.v1.x*(m.v2.y*m.v3.z - m.v2.z*m.v3.y) +
        m.v2.x*(m.v1.z*m.v3.y - m.v1.y*m.v3.z) +
        m.v3.x*(m.v1.y*m.v2.z - m.v1.z*m.v2.y));

    out.v2.x = (m.v2.z*m.v3.x - m.v2.x*m.v3.z) / (
        m.v1.x*(m.v2.y*m.v3.z - m.v2.z*m.v3.y) +
        m.v2.x*(m.v1.z*m.v3.y - m.v1.y*m.v3.z) +
        m.v3.x*(m.v1.y*m.v2.z - m.v1.z*m.v2.y));

    out.v3.x = (m.v2.x*m.v3.y - m.v2.y*m.v3.x) / (
        m.v1.x*(m.v2.y*m.v3.z - m.v2.z*m.v3.y) +
        m.v2.x*(m.v1.z*m.v3.y - m.v1.y*m.v3.z) +
        m.v3.x*(m.v1.y*m.v2.z - m.v1.z*m.v2.y));

    out.v1.y = (m.v1.z*m.v3.y - m.v1.y*m.v3.z) / (
        m.v1.x*(m.v2.y*m.v3.z - m.v2.z*m.v3.y) +
        m.v2.x*(m.v1.z*m.v3.y - m.v1.y*m.v3.z) +
        m.v3.x*(m.v1.y*m.v2.z - m.v1.z*m.v2.y));

    out.v2.y = (m.v1.x*m.v3.z - m.v1.z*m.v3.x) / (
        m.v1.x*(m.v2.y*m.v3.z - m.v2.z*m.v3.y) +
        m.v2.x*(m.v1.z*m.v3.y - m.v1.y*m.v3.z) +
        m.v3.x*(m.v1.y*m.v2.z - m.v1.z*m.v2.y));

    out.v3.y = (m.v1.y*m.v3.x - m.v1.x*m.v3.y) / (
        m.v1.x*(m.v2.y*m.v3.z - m.v2.z*m.v3.y) +
        m.v2.x*(m.v1.z*m.v3.y - m.v1.y*m.v3.z) +
        m.v3.x*(m.v1.y*m.v2.z - m.v1.z*m.v2.y));

    out.v1.z = (m.v1.y*m.v2.z - m.v1.z*m.v2.y) / (
        m.v1.x*(m.v2.y*m.v3.z - m.v2.z*m.v3.y) +
        m.v2.x*(m.v1.z*m.v3.y - m.v1.y*m.v3.z) +
        m.v3.x*(m.v1.y*m.v2.z - m.v1.z*m.v2.y));

    out.v2.z = (m.v1.z*m.v2.x - m.v1.x*m.v2.z) / (
        m.v1.x*(m.v2.y*m.v3.z - m.v2.z*m.v3.y) +
        m.v2.x*(m.v1.z*m.v3.y - m.v1.y*m.v3.z) +
        m.v3.x*(m.v1.y*m.v2.z - m.v1.z*m.v2.y));

    out.v3.z = (m.v1.x*m.v2.y - m.v1.y*m.v2.x) / (
        m.v1.x*(m.v2.y*m.v3.z - m.v2.z*m.v3.y) +
        m.v2.x*(m.v1.z*m.v3.y - m.v1.y*m.v3.z) +
        m.v3.x*(m.v1.y*m.v2.z - m.v1.z*m.v2.y));

    return out;
  }

  mat4 inverse(mat4 m)
  {
    static mat4 out;

    out.v1.x = (
        m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
        m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
        m.v4.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z)
        ) / (
          m.v1.x * (
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            m.v4.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z)
            ) +
          m.v2.x * (
            -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
            m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
            m.v4.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z)
            ) + m.v3.x * (
              m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
              m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
              m.v4.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
              ) + m.v4.x * (
                -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                m.v3.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
                )
          );

    out.v2.x = (
        -m.v2.x*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
        m.v3.x*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
        m.v4.x*(m.v2.z*m.v3.w - m.v2.w*m.v3.z)
        ) / (
          m.v1.x * (
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            m.v4.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z)
            ) + m.v2.x * (
              -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
              m.v4.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z)
              ) + m.v3.x * (
                m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
                m.v4.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
                ) + m.v4.x * (
                  -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                  m.v3.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
                  )
          );

    out.v3.x = (
        m.v2.x*(m.v3.y*m.v4.w - m.v3.w*m.v4.y) -
        m.v3.x*(m.v2.y*m.v4.w - m.v2.w*m.v4.y) +
        m.v4.x*(m.v2.y*m.v3.w - m.v2.w*m.v3.y)
        ) / (
          m.v1.x * (
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            m.v4.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z)
            ) + m.v2.x * (
              -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
              m.v4.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z)
              ) + m.v3.x * (
                m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
                m.v4.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
                ) + m.v3.x * (
                  -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                  m.v3.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
                  )
          );

    out .v4.x = (
        -m.v2.x*(m.v3.y*m.v4.z - m.v3.z*m.v4.y) +
        m.v3.x*(m.v2.y*m.v4.z - m.v2.z*m.v4.y) -
        m.v4.x*(m.v2.y*m.v3.z - m.v2.z*m.v3.y)
        ) / (
          m.v1.x * (
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            (m.v2.z*m.v3.w - m.v2.w*m.v3.z)*m.v4.y
            ) + m.v2.x * (
              -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
              m.v4.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z)
              ) + m.v3.x * (
                m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
                m.v4.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
                ) + m.v4.x * (
                  -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                  m.v3.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
                  )
          );

    out.v1.y = (
        -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
        m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
        m.v4.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z)
        ) / (
          m.v1.x * (
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            m.v4.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z)
            ) + m.v2.x * (
              -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
              m.v4.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z)
              ) + m.v3.x * (
                m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
                m.v4.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
                ) + m.v4.x * (
                  -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                  m.v3.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
                  )
          );

    out.v2.y = (
        m.v1.x*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
        m.v3.x*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
        m.v4.x*(m.v1.z*m.v3.w - m.v1.w*m.v3.z)
        ) / (
          m.v1.x * (
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            m.v4.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z)
            ) + m.v2.x * (
              -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
              m.v4.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z)
              ) + m.v3.x * (
                m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
                m.v4.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
                ) + m.v4.x * (
                  -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                  m.v3.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
                  )
          );

    out.v3.y = (
        -m.v1.x*(m.v3.y*m.v4.w - m.v3.w*m.v4.y) +
        m.v3.x*(m.v1.y*m.v4.w - m.v1.w*m.v4.y) -
        m.v4.x*(m.v1.y*m.v3.w - m.v1.w*m.v3.y)
        ) / (
          m.v1.x * (
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            m.v4.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z)
            ) + m.v2.x * (
              -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
              m.v4.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z)
              ) + m.v3.x * (
                m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
                m.v4.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
                ) + m.v4.x * (
                  -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                  m.v3.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
                  )
          );

    out.v4.y = (
        m.v1.x*(m.v3.y*m.v4.z - m.v3.z*m.v4.y) -
        m.v3.x*(m.v1.y*m.v4.z - m.v1.z*m.v4.y) +
        (m.v1.y*m.v3.z - m.v1.z*m.v3.y)*m.v4.x
        ) / (
          m.v1.x * (
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            (m.v2.z*m.v3.w - m.v2.w*m.v3.z)*m.v4.y
            ) + m.v2.x * (
              -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
              (m.v1.z*m.v3.w - m.v1.w*m.v3.z)*m.v4.y
              ) + m.v3.x * (
                m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
                (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v4.y
                ) + m.v4.x * (
                  -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                  (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v3.y
                  )
          );

    out.v1.z = (
        m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
        m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
        (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v4.y
        ) / (
          m.v1.x * (
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            (m.v2.z*m.v3.w - m.v2.w*m.v3.z)*m.v4.y
            ) + m.v2.x * (
              -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
              (m.v1.z*m.v3.w - m.v1.w*m.v3.z)*m.v4.y
              ) + m.v3.x * (
                m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
                (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v4.y
                ) + m.v4.x * (
                  -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                  (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v3.y
                  )
          );

    out.v2.z = (
        -m.v1.x*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
        m.v2.x*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
        (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v4.x
        ) / (
          m.v1.x*(
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            (m.v2.z*m.v3.w - m.v2.w*m.v3.z)*m.v4.y
            ) + m.v2.x*(
              -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
              (m.v1.z*m.v3.w - m.v1.w*m.v3.z)*m.v4.y
              ) + m.v3.x*(
                m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
                (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v4.y
                ) + m.v4.x*(
                  -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                  (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v3.y
                  )
          );

    out.v3.z = (
        m.v1.x*(m.v2.y*m.v4.w - m.v2.w*m.v4.y) -
        m.v2.x*(m.v1.y*m.v4.w - m.v1.w*m.v4.y) +
        (m.v1.y*m.v2.w - m.v1.w*m.v2.y)*m.v4.x
        ) / (
          m.v1.x*(
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            (m.v2.z*m.v3.w - m.v2.w*m.v3.z)*m.v4.y
            ) + m.v2.x*(
              -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
              (m.v1.z*m.v3.w - m.v1.w*m.v3.z)*m.v4.y
              ) + m.v3.x*(
                m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
                (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v4.y
                ) + m.v4.x*(
                  -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                  (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v3.y
                  )
          );

    out.v4.z = (
        -m.v1.x*(m.v2.y*m.v4.z - m.v2.z*m.v4.y) +
        m.v2.x*(m.v1.y*m.v4.z - m.v1.z*m.v4.y) -
        (m.v1.y*m.v2.z - m.v1.z*m.v2.y)*m.v4.x
        ) / (
          m.v1.x*(
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            (m.v2.z*m.v3.w - m.v2.w*m.v3.z)*m.v4.y
            ) + m.v2.x*(
              -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
              (m.v1.z*m.v3.w - m.v1.w*m.v3.z)*m.v4.y
              ) + m.v3.x*(
                m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
                (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v4.y
                ) + m.v4.x*(
                  -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                  (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v3.y
                  )
          );

    out.v1.w = (
        -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
        m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
        (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v3.y
        ) / (
          m.v1.x*(
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            (m.v2.z*m.v3.w - m.v2.w*m.v3.z)*m.v4.y
            ) + m.v2.x*(
              -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
              (m.v1.z*m.v3.w - m.v1.w*m.v3.z)*m.v4.y
              ) + m.v3.x*(
                m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
                (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v4.y
                ) + (
                  -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                  (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v3.y
                  )*m.v4.x
          );

    out.v2.w = (
        m.v1.x*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) -
        m.v2.x*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) +
        (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v3.x
        ) / (
          m.v1.x*(
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            (m.v2.z*m.v3.w - m.v2.w*m.v3.z)*m.v4.y
            ) + m.v2.x*(
              -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
              (m.v1.z*m.v3.w - m.v1.w*m.v3.z)*m.v4.y
              ) + m.v3.x*(
                m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
                (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v4.y
                ) + m.v4.x*(
                  -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                  (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v3.y
                  )
          );

    out.v3.w = (
        -m.v1.x*(m.v2.y*m.v3.w - m.v2.w*m.v3.y) +
        m.v2.x*(m.v1.y*m.v3.w - m.v1.w*m.v3.y) -
        (m.v1.y*m.v2.w - m.v1.w*m.v2.y)*m.v3.x
        ) / (
          m.v1.x*(
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            (m.v2.z*m.v3.w - m.v2.w*m.v3.z)*m.v4.y
            ) + m.v2.x*(
              -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
              (m.v1.z*m.v3.w - m.v1.w*m.v3.z)*m.v4.y
              ) + m.v3.x*(
                m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
                (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v4.y
                ) + m.v4.x*(
                  -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                  (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v3.y
                  )
          );

    out.v4.w = (
        m.v1.x*(m.v2.y*m.v3.z-m.v2.z*m.v3.y) -
        m.v2.x*(m.v1.y*m.v3.z-m.v1.z*m.v3.y) +
        (m.v1.y*m.v2.z-m.v1.z*m.v2.y)*m.v3.x
        ) / (
          m.v1.x*(
            m.v2.y*(m.v3.z*m.v4.w-m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w-m.v2.w*m.v4.z) +
            (m.v2.z*m.v3.w-m.v2.w*m.v3.z)*m.v4.y
            ) + m.v2.x*(
              -m.v1.y*(m.v3.z*m.v4.w-m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w-m.v1.w*m.v4.z) -
              (m.v1.z*m.v3.w-m.v1.w*m.v3.z)*m.v4.y
              ) + m.v3.x*(
                m.v1.y*(m.v2.z*m.v4.w-m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w-m.v1.w*m.v4.z) +
                (m.v1.z*m.v2.w-m.v1.w*m.v2.z)*m.v4.y
                ) + m.v4.x*(
                  -m.v1.y*(m.v2.z*m.v3.w-m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w-m.v1.w*m.v3.z) -
                  (m.v1.z*m.v2.w-m.v1.w*m.v2.z)*m.v3.y
                  )
          );

    return out;
  }

  dmat2 inverse(dmat2 m)
  {
    static dmat2 out;

    out.v1.x =  m.v2.y/(m.v1.x*m.v2.y - m.v1.y*m.v2.x);
    out.v1.y = -m.v1.y/(m.v1.x*m.v2.y - m.v1.y*m.v2.x);
    out.v2.x = -m.v2.x/(m.v1.x*m.v2.y - m.v1.y*m.v2.x);
    out.v2.y =  m.v1.x/(m.v1.x*m.v2.y - m.v1.y*m.v2.x);

    return out;
  }

  dmat3 inverse(dmat3 m)
  {
    static dmat3 out;

    out.v1.x = (m.v2.y*m.v3.z - m.v2.z*m.v3.y) / (
        m.v1.x*(m.v2.y*m.v3.z - m.v2.z*m.v3.y) +
        m.v2.x*(m.v1.z*m.v3.y - m.v1.y*m.v3.z) +
        m.v3.x*(m.v1.y*m.v2.z - m.v1.z*m.v2.y));

    out.v2.x = (m.v2.z*m.v3.x - m.v2.x*m.v3.z) / (
        m.v1.x*(m.v2.y*m.v3.z - m.v2.z*m.v3.y) +
        m.v2.x*(m.v1.z*m.v3.y - m.v1.y*m.v3.z) +
        m.v3.x*(m.v1.y*m.v2.z - m.v1.z*m.v2.y));

    out.v3.x = (m.v2.x*m.v3.y - m.v2.y*m.v3.x) / (
        m.v1.x*(m.v2.y*m.v3.z - m.v2.z*m.v3.y) +
        m.v2.x*(m.v1.z*m.v3.y - m.v1.y*m.v3.z) +
        m.v3.x*(m.v1.y*m.v2.z - m.v1.z*m.v2.y));

    out.v1.y = (m.v1.z*m.v3.y - m.v1.y*m.v3.z) / (
        m.v1.x*(m.v2.y*m.v3.z - m.v2.z*m.v3.y) +
        m.v2.x*(m.v1.z*m.v3.y - m.v1.y*m.v3.z) +
        m.v3.x*(m.v1.y*m.v2.z - m.v1.z*m.v2.y));

    out.v2.y = (m.v1.x*m.v3.z - m.v1.z*m.v3.x) / (
        m.v1.x*(m.v2.y*m.v3.z - m.v2.z*m.v3.y) +
        m.v2.x*(m.v1.z*m.v3.y - m.v1.y*m.v3.z) +
        m.v3.x*(m.v1.y*m.v2.z - m.v1.z*m.v2.y));

    out.v3.y = (m.v1.y*m.v3.x - m.v1.x*m.v3.y) / (
        m.v1.x*(m.v2.y*m.v3.z - m.v2.z*m.v3.y) +
        m.v2.x*(m.v1.z*m.v3.y - m.v1.y*m.v3.z) +
        m.v3.x*(m.v1.y*m.v2.z - m.v1.z*m.v2.y));

    out.v1.z = (m.v1.y*m.v2.z - m.v1.z*m.v2.y) / (
        m.v1.x*(m.v2.y*m.v3.z - m.v2.z*m.v3.y) +
        m.v2.x*(m.v1.z*m.v3.y - m.v1.y*m.v3.z) +
        m.v3.x*(m.v1.y*m.v2.z - m.v1.z*m.v2.y));

    out.v2.z = (m.v1.z*m.v2.x - m.v1.x*m.v2.z) / (
        m.v1.x*(m.v2.y*m.v3.z - m.v2.z*m.v3.y) +
        m.v2.x*(m.v1.z*m.v3.y - m.v1.y*m.v3.z) +
        m.v3.x*(m.v1.y*m.v2.z - m.v1.z*m.v2.y));

    out.v3.z = (m.v1.x*m.v2.y - m.v1.y*m.v2.x) / (
        m.v1.x*(m.v2.y*m.v3.z - m.v2.z*m.v3.y) +
        m.v2.x*(m.v1.z*m.v3.y - m.v1.y*m.v3.z) +
        m.v3.x*(m.v1.y*m.v2.z - m.v1.z*m.v2.y));

    return out;
  }

  dmat4 inverse(dmat4 m)
  {
    static dmat4 out;

    out.v1.x = (
        m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
        m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
        m.v4.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z)
        ) / (
          m.v1.x * (
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            m.v4.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z)
            ) +
          m.v2.x * (
            -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
            m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
            m.v4.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z)
            ) + m.v3.x * (
              m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
              m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
              m.v4.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
              ) + m.v4.x * (
                -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                m.v3.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
                )
          );

    out.v2.x = (
        -m.v2.x*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
        m.v3.x*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
        m.v4.x*(m.v2.z*m.v3.w - m.v2.w*m.v3.z)
        ) / (
          m.v1.x * (
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            m.v4.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z)
            ) + m.v2.x * (
              -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
              m.v4.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z)
              ) + m.v3.x * (
                m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
                m.v4.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
                ) + m.v4.x * (
                  -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                  m.v3.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
                  )
          );

    out.v3.x = (
        m.v2.x*(m.v3.y*m.v4.w - m.v3.w*m.v4.y) -
        m.v3.x*(m.v2.y*m.v4.w - m.v2.w*m.v4.y) +
        m.v4.x*(m.v2.y*m.v3.w - m.v2.w*m.v3.y)
        ) / (
          m.v1.x * (
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            m.v4.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z)
            ) + m.v2.x * (
              -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
              m.v4.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z)
              ) + m.v3.x * (
                m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
                m.v4.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
                ) + m.v3.x * (
                  -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                  m.v3.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
                  )
          );

    out .v4.x = (
        -m.v2.x*(m.v3.y*m.v4.z - m.v3.z*m.v4.y) +
        m.v3.x*(m.v2.y*m.v4.z - m.v2.z*m.v4.y) -
        m.v4.x*(m.v2.y*m.v3.z - m.v2.z*m.v3.y)
        ) / (
          m.v1.x * (
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            (m.v2.z*m.v3.w - m.v2.w*m.v3.z)*m.v4.y
            ) + m.v2.x * (
              -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
              m.v4.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z)
              ) + m.v3.x * (
                m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
                m.v4.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
                ) + m.v4.x * (
                  -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                  m.v3.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
                  )
          );

    out.v1.y = (
        -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
        m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
        m.v4.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z)
        ) / (
          m.v1.x * (
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            m.v4.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z)
            ) + m.v2.x * (
              -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
              m.v4.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z)
              ) + m.v3.x * (
                m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
                m.v4.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
                ) + m.v4.x * (
                  -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                  m.v3.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
                  )
          );

    out.v2.y = (
        m.v1.x*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
        m.v3.x*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
        m.v4.x*(m.v1.z*m.v3.w - m.v1.w*m.v3.z)
        ) / (
          m.v1.x * (
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            m.v4.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z)
            ) + m.v2.x * (
              -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
              m.v4.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z)
              ) + m.v3.x * (
                m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
                m.v4.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
                ) + m.v4.x * (
                  -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                  m.v3.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
                  )
          );

    out.v3.y = (
        -m.v1.x*(m.v3.y*m.v4.w - m.v3.w*m.v4.y) +
        m.v3.x*(m.v1.y*m.v4.w - m.v1.w*m.v4.y) -
        m.v4.x*(m.v1.y*m.v3.w - m.v1.w*m.v3.y)
        ) / (
          m.v1.x * (
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            m.v4.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z)
            ) + m.v2.x * (
              -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
              m.v4.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z)
              ) + m.v3.x * (
                m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
                m.v4.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
                ) + m.v4.x * (
                  -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                  m.v3.y*(m.v1.z*m.v2.w - m.v1.w*m.v2.z)
                  )
          );

    out.v4.y = (
        m.v1.x*(m.v3.y*m.v4.z - m.v3.z*m.v4.y) -
        m.v3.x*(m.v1.y*m.v4.z - m.v1.z*m.v4.y) +
        (m.v1.y*m.v3.z - m.v1.z*m.v3.y)*m.v4.x
        ) / (
          m.v1.x * (
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            (m.v2.z*m.v3.w - m.v2.w*m.v3.z)*m.v4.y
            ) + m.v2.x * (
              -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
              (m.v1.z*m.v3.w - m.v1.w*m.v3.z)*m.v4.y
              ) + m.v3.x * (
                m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
                (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v4.y
                ) + m.v4.x * (
                  -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                  (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v3.y
                  )
          );

    out.v1.z = (
        m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
        m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
        (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v4.y
        ) / (
          m.v1.x * (
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            (m.v2.z*m.v3.w - m.v2.w*m.v3.z)*m.v4.y
            ) + m.v2.x * (
              -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
              (m.v1.z*m.v3.w - m.v1.w*m.v3.z)*m.v4.y
              ) + m.v3.x * (
                m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
                (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v4.y
                ) + m.v4.x * (
                  -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                  (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v3.y
                  )
          );

    out.v2.z = (
        -m.v1.x*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
        m.v2.x*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
        (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v4.x
        ) / (
          m.v1.x*(
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            (m.v2.z*m.v3.w - m.v2.w*m.v3.z)*m.v4.y
            ) + m.v2.x*(
              -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
              (m.v1.z*m.v3.w - m.v1.w*m.v3.z)*m.v4.y
              ) + m.v3.x*(
                m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
                (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v4.y
                ) + m.v4.x*(
                  -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                  (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v3.y
                  )
          );

    out.v3.z = (
        m.v1.x*(m.v2.y*m.v4.w - m.v2.w*m.v4.y) -
        m.v2.x*(m.v1.y*m.v4.w - m.v1.w*m.v4.y) +
        (m.v1.y*m.v2.w - m.v1.w*m.v2.y)*m.v4.x
        ) / (
          m.v1.x*(
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            (m.v2.z*m.v3.w - m.v2.w*m.v3.z)*m.v4.y
            ) + m.v2.x*(
              -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
              (m.v1.z*m.v3.w - m.v1.w*m.v3.z)*m.v4.y
              ) + m.v3.x*(
                m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
                (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v4.y
                ) + m.v4.x*(
                  -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                  (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v3.y
                  )
          );

    out.v4.z = (
        -m.v1.x*(m.v2.y*m.v4.z - m.v2.z*m.v4.y) +
        m.v2.x*(m.v1.y*m.v4.z - m.v1.z*m.v4.y) -
        (m.v1.y*m.v2.z - m.v1.z*m.v2.y)*m.v4.x
        ) / (
          m.v1.x*(
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            (m.v2.z*m.v3.w - m.v2.w*m.v3.z)*m.v4.y
            ) + m.v2.x*(
              -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
              (m.v1.z*m.v3.w - m.v1.w*m.v3.z)*m.v4.y
              ) + m.v3.x*(
                m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
                (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v4.y
                ) + m.v4.x*(
                  -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                  (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v3.y
                  )
          );

    out.v1.w = (
        -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
        m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
        (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v3.y
        ) / (
          m.v1.x*(
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            (m.v2.z*m.v3.w - m.v2.w*m.v3.z)*m.v4.y
            ) + m.v2.x*(
              -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
              (m.v1.z*m.v3.w - m.v1.w*m.v3.z)*m.v4.y
              ) + m.v3.x*(
                m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
                (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v4.y
                ) + (
                  -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                  (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v3.y
                  )*m.v4.x
          );

    out.v2.w = (
        m.v1.x*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) -
        m.v2.x*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) +
        (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v3.x
        ) / (
          m.v1.x*(
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            (m.v2.z*m.v3.w - m.v2.w*m.v3.z)*m.v4.y
            ) + m.v2.x*(
              -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
              (m.v1.z*m.v3.w - m.v1.w*m.v3.z)*m.v4.y
              ) + m.v3.x*(
                m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
                (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v4.y
                ) + m.v4.x*(
                  -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                  (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v3.y
                  )
          );

    out.v3.w = (
        -m.v1.x*(m.v2.y*m.v3.w - m.v2.w*m.v3.y) +
        m.v2.x*(m.v1.y*m.v3.w - m.v1.w*m.v3.y) -
        (m.v1.y*m.v2.w - m.v1.w*m.v2.y)*m.v3.x
        ) / (
          m.v1.x*(
            m.v2.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) +
            (m.v2.z*m.v3.w - m.v2.w*m.v3.z)*m.v4.y
            ) + m.v2.x*(
              -m.v1.y*(m.v3.z*m.v4.w - m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) -
              (m.v1.z*m.v3.w - m.v1.w*m.v3.z)*m.v4.y
              ) + m.v3.x*(
                m.v1.y*(m.v2.z*m.v4.w - m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w - m.v1.w*m.v4.z) +
                (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v4.y
                ) + m.v4.x*(
                  -m.v1.y*(m.v2.z*m.v3.w - m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w - m.v1.w*m.v3.z) -
                  (m.v1.z*m.v2.w - m.v1.w*m.v2.z)*m.v3.y
                  )
          );

    out.v4.w = (
        m.v1.x*(m.v2.y*m.v3.z-m.v2.z*m.v3.y) -
        m.v2.x*(m.v1.y*m.v3.z-m.v1.z*m.v3.y) +
        (m.v1.y*m.v2.z-m.v1.z*m.v2.y)*m.v3.x
        ) / (
          m.v1.x*(
            m.v2.y*(m.v3.z*m.v4.w-m.v3.w*m.v4.z) -
            m.v3.y*(m.v2.z*m.v4.w-m.v2.w*m.v4.z) +
            (m.v2.z*m.v3.w-m.v2.w*m.v3.z)*m.v4.y
            ) + m.v2.x*(
              -m.v1.y*(m.v3.z*m.v4.w-m.v3.w*m.v4.z) +
              m.v3.y*(m.v1.z*m.v4.w-m.v1.w*m.v4.z) -
              (m.v1.z*m.v3.w-m.v1.w*m.v3.z)*m.v4.y
              ) + m.v3.x*(
                m.v1.y*(m.v2.z*m.v4.w-m.v2.w*m.v4.z) -
                m.v2.y*(m.v1.z*m.v4.w-m.v1.w*m.v4.z) +
                (m.v1.z*m.v2.w-m.v1.w*m.v2.z)*m.v4.y
                ) + m.v4.x*(
                  -m.v1.y*(m.v2.z*m.v3.w-m.v2.w*m.v3.z) +
                  m.v2.y*(m.v1.z*m.v3.w-m.v1.w*m.v3.z) -
                  (m.v1.z*m.v2.w-m.v1.w*m.v2.z)*m.v3.y
                  )
          );

    return out;
  }

  mat3 translate(mat3 m, float x, float y)
  {
    static mat3 out(1.0f);

    out.v3.x = x;
    out.v3.y = y;

    out = out * m;

    return out;
  }

  mat3 translate(mat3 m, vec2 v)
  {
    return translate(m, v.x, v.y);
  }

  mat4 translate(mat4 m, float x, float y, float z)
  {
    static mat4 out(1.0f);

    out.v4.x = x;
    out.v4.y = y;
    out.v4.z = z;

    out = out * m;

    return out;
  }

  mat4 translate(mat4 m, vec3 v)
  {
    return translate(m, v.x, v.y, v.z);
  }

  dmat3 translate(dmat3 m, double x, double y)
  {
    static dmat3 out(1.0f);

    out.v3.x = x;
    out.v3.y = y;

    out = out * m;

    return out;
  }

  dmat3 translate(dmat3 m, dvec2 v)
  {
    return translate(m, v.x, v.y);
  }

  dmat4 translate(dmat4 m, double x, double y, double z)
  {
    static dmat4 out(1.0f);

    out.v4.x = x;
    out.v4.y = y;
    out.v4.z = z;

    out = out * m;

    return out;
  }

  dmat4 translate(dmat4 m, dvec3 v)
  {
    return translate(m, v.x, v.y, v.z);
  }

  mat2 scaleX(mat2 m, float k)
  {
    static mat2 out(1.0f);

    out.v1.x = k;

    out = out * m;

    return out;
  }

  mat3 scaleX(mat3 m, float k)
  {
    static mat3 out(1.0f);

    out.v1.x = k;

    out = out * m;

    return out;
  }

  mat4 scaleX(mat4 m, float k)
  {
    static mat4 out(1.0f);

    out.v1.x = k;

    out = out * m;

    return out;
  }

  mat3 scaleY(mat3 m, float k)
  {
    static mat3 out(1.0f);

    out.v2.y = k;

    out = out * m;

    return out;
  }

  mat4 scaleY(mat4 m, float k)
  {
    static mat4 out(1.0f);

    out.v2.y = k;

    out = out * m;

    return out;
  }

  mat4 scaleZ(mat4 m, float k)
  {
    static mat4 out(1.0f);

    out.v3.z = k;

    out = out * m;

    return out;
  }

  dmat3 scaleX(dmat3 m, double k)
  {
    static dmat3 out(1.0f);

    out.v1.x = k;

    out = out * m;

    return out;
  }

  dmat4 scaleX(dmat4 m, double k)
  {
    static dmat4 out(1.0f);

    out.v1.x = k;

    out = out * m;

    return out;
  }

  dmat3 scaleY(dmat3 m, double k)
  {
    static dmat3 out(1.0f);

    out.v2.y = k;

    out = out * m;

    return out;
  }

  dmat4 scaleY(dmat4 m, double k)
  {
    static dmat4 out(1.0f);

    out.v2.y = k;

    out = out * m;

    return out;
  }

  dmat4 scaleZ(dmat4 m, double k)
  {
    static dmat4 out(1.0f);

    out.v3.z = k;

    out = out * m;

    return out;
  }

  mat3 scale(mat3 m, float k)
  {
    static mat3 out(1.0f);

    out.v1.x = k;
    out.v2.y = k;
    out.v3.z = k;

    out = out * m;

    return out;
  }

  mat4 scale(mat4 m, float k)
  {
    static mat4 out(1.0f);

    out.v1.x = k;
    out.v2.y = k;
    out.v3.z = k;

    out = out * m;

    return out;
  }

  dmat3 scale(dmat3 m, double k)
  {
    static dmat3 out(1.0f);

    out.v1.x = k;
    out.v2.y = k;
    out.v3.z = k;

    out = out * m;

    return out;
  }

  dmat4 scale(dmat4 m, double k)
  {
    static dmat4 out(1.0f);

    out.v1.x = k;
    out.v2.y = k;
    out.v3.z = k;

    out = out * m;

    return out;
  }

  mat4 rotateX(mat4 m, float angle)
  {
    static mat4 out(1.0f);

    float s = sinf(angle);
    float c = cosf(angle);

    out.v2.y = c;
    out.v2.z = s;
    out.v3.y = -s;
    out.v3.z = c;

    out = out * m;

    return out;
  }

  mat4 rotateY(mat4 m, float angle)
  {
    static mat4 out(1.0f);

    float s = sinf(angle);
    float c = cosf(angle);

    out.v1.x = c;
    out.v1.z = -s;
    out.v3.x = s;
    out.v3.z = c;

    out = out * m;

    return out;
  }

  mat4 rotateZ(mat4 m, float angle)
  {
    static mat4 out(1.0f);

    float s = sinf(angle);
    float c = cosf(angle);

    out.v1.x = c;
    out.v1.y = s;
    out.v2.x = -s;
    out.v2.y = c;

    out = out * m;

    return out;
  }

  dmat4 rotateX(dmat4 m, double angle)
  {
    static dmat4 out(1.0f);

    double s = sin(angle);
    double c = cos(angle);

    out.v2.y = c;
    out.v2.z = s;
    out.v3.y = -s;
    out.v3.z = c;

    out = out * m;

    return out;
  }

  dmat4 rotateY(dmat4 m, double angle)
  {
    static dmat4 out(1.0f);

    double s = sin(angle);
    double c = cos(angle);

    out.v1.x = c;
    out.v1.z = -s;
    out.v3.x = s;
    out.v3.z = c;

    out = out * m;

    return out;
  }

  dmat4 rotateZ(dmat4 m, double angle)
  {
    static dmat4 out(1.0f);

    double s = sin(angle);
    double c = cos(angle);

    out.v1.x = c;
    out.v1.y = s;
    out.v2.x = -s;
    out.v2.y = c;

    out = out * m;

    return out;
  }

  mat3 rotate(mat3 m, float angle)
  {
    static mat3 out(1.0f);

    double s = sinf(angle);
    double c = cosf(angle);

    out.v1.x = c;
    out.v1.y = s;
    out.v2.x = -s;
    out.v2.y = c;

    out = out * m;

    return out;
  }

  mat4 rotate(mat4 m, float angle, vec3 v)
  {
    return rotate(m, angle, v.x, v.y, v.z);
  }

  mat4 rotate(mat4 m, float angle, float x, float y, float z)
  {
    static mat4 out;

    vec3 a(x, y, z);
    a = normalize(a);

    float sv = sinf(angle);
    float cv = cosf(angle);
    vec3 mv = a * (1 - cv);

    out.v1.x = a.x * mv.x + cv;
    out.v1.y = a.x * mv.y + a.z * sv;
    out.v1.z = a.x * mv.z - a.y * sv;

    out.v2.x = a.y * mv.x - a.z * sv;
    out.v2.y = a.y * mv.y + cv;
    out.v2.z = a.y * mv.z + a.x * sv;

    out.v3.x = a.z * mv.x + a.y * sv;
    out.v3.y = a.z * mv.y - a.x * sv;
    out.v3.z = a.z * mv.z + cv;

    out.v4.w = 1.0f;

    out = out * m;

    return out;
  }

  dmat3 rotate(dmat3 m, double angle)
  {
    static dmat3 out(1.0f);

    double s = sin(angle);
    double c = cos(angle);

    out.v1.x = c;
    out.v1.y = s;
    out.v2.x = -s;
    out.v2.y = c;

    out = out * m;

    return out;
  }

  dmat4 rotate(dmat4 m, double angle, dvec3 v)
  {
    return rotate(m, angle, v.x, v.y, v.z);
  }

  dmat4 rotate(dmat4 m, double angle, double x, double y, double z)
  {
    static dmat4 out;

    dvec3 a(x, y, z);
    a = normalize(a);

    double sv = sin(angle);
    double cv = cos(angle);
    dvec3 mv = a * (1 - cv);

    out.v1.x = a.x * mv.x + cv;
    out.v1.y = a.x * mv.y + a.z * sv;
    out.v1.z = a.x * mv.z - a.y * sv;

    out.v2.x = a.y * mv.x - a.z * sv;
    out.v2.y = a.y * mv.y + cv;
    out.v2.z = a.y * mv.z + a.x * sv;

    out.v3.x = a.z * mv.x + a.y * sv;
    out.v3.y = a.z * mv.y - a.x * sv;
    out.v3.z = a.z * mv.z + cv;

    out.v4.w = 1.0f;

    out = out * m;

    return out;
  }

  mat4 orthof(float left, float right, float bottom, float top, float near, float far)
  {
    static mat4 out;

    out.v1.x = 2.0f / (right - left);
    out.v1.y = 0.0f;
    out.v1.z = 0.0f;
    out.v1.w = 0.0f;

    out.v2.x = 0.0f;
    out.v2.y = 2.0f / (top - bottom);
    out.v2.z = 0.0f;
    out.v2.w = 0.0f;

    out.v3.x = 0.0f;
    out.v3.y = 0.0f;
    out.v3.z = -2.0f / (far - near);
    out.v3.w = 0.0f;

    out.v4.x = -(right + left) / (right - left);
    out.v4.y = -(top + bottom) / (top - bottom);
    out.v4.z = -(far + near) / (far - near);
    out.v4.w = 1.0f;

    return out;
  }

  dmat4 orthod(double left, double right, double bottom, double top, double near, double far)
  {
    static dmat4 out;

    out.v1.x = 2.0f / (right - left);
    out.v1.y = 0.0f;
    out.v1.z = 0.0f;
    out.v1.w = 0.0f;

    out.v2.x = 0.0f;
    out.v2.y = 2.0f / (top - bottom);
    out.v2.z = 0.0f;
    out.v2.w = 0.0f;

    out.v3.x = 0.0f;
    out.v3.y = 0.0f;
    out.v3.z = -2.0f / (far - near);
    out.v3.w = 0.0f;

    out.v4.x = -(right + left) / (right - left);
    out.v4.y = -(top + bottom) / (top - bottom);
    out.v4.z = -(far + near) / (far - near);
    out.v4.w = 1.0f;

    return out;
  }

  mat4 frustumf(float left, float right, float bottom, float top, float near, float far)
  {
    static mat4 out;

    out.v1.x = 2.0f / (right - left);
    out.v1.y = 0.0f;
    out.v1.z = 0.0f;
    out.v1.w = 0.0f;

    out.v2.x = 0.0f;
    out.v2.y = 2.0f / (top - bottom);
    out.v2.z = 0.0f;
    out.v2.w = 0.0f;

    out.v3.x = (right + left) / (right - left);
    out.v3.y = (top + bottom) / (top - bottom);
    out.v3.z = -(far + near) / (far - near);
    out.v3.w = -1.0f;

    out.v4.x = 0.0f;
    out.v4.y = 0.0f;
    out.v4.z = -2.0f * (far * near) / (far - near);
    out.v4.w = 0.0f;

    return out;
  }

  dmat4 frustumd(double left, double right, double bottom, double top, double near, double far)
  {
    static dmat4 out;

    out.v1.x = 2.0f / (right - left);
    out.v1.y = 0.0f;
    out.v1.z = 0.0f;
    out.v1.w = 0.0f;

    out.v2.x = 0.0f;
    out.v2.y = 2.0f / (top - bottom);
    out.v2.z = 0.0f;
    out.v2.w = 0.0f;

    out.v3.x = (right + left) / (right - left);
    out.v3.y = (top + bottom) / (top - bottom);
    out.v3.z = -(far + near) / (far - near);
    out.v3.w = -1.0f;

    out.v4.x = 0.0f;
    out.v4.y = 0.0f;
    out.v4.z = -2.0f * (far * near) / (far - near);
    out.v4.w = 0.0f;

    return out;
  }

  mat4 perspectivef(float fov, float aspect, float near, float far)
  {
    static mat4 out;

    float a = 1.0f / tan(fov / 2.0f);

    out.v1.x = a / aspect;
    out.v1.y = 0.0f;
    out.v1.z = 0.0f;
    out.v1.w = 0.0f;

    out.v2.x = 0.0f;
    out.v2.y = a;
    out.v2.z = 0.0f;
    out.v2.w = 0.0f;

    out.v3.x = 0.0f;
    out.v3.y = 0.0f;
    out.v3.z = -(far + near) / (far - near);
    out.v3.w = -1.0f;

    out.v4.x = 0.0f;
    out.v4.y = 0.0f;
    out.v4.z = -(2.0f * far * near) / (far - near);
    out.v4.w = 0.0f;

    return out;
  }

  dmat4 perspectived(double fov, double aspect, double near, double far)
  {
    static dmat4 out;

    double a = 1.0f / tan(fov / 2.0f);

    out.v1.x = a / aspect;
    out.v1.y = 0.0f;
    out.v1.z = 0.0f;
    out.v1.w = 0.0f;

    out.v2.x = 0.0f;
    out.v2.y = a;
    out.v2.z = 0.0f;
    out.v2.w = 0.0f;

    out.v3.x = 0.0f;
    out.v3.y = 0.0f;
    out.v3.z = -(far + near) / (far - near);
    out.v3.w = -1.0f;

    out.v4.x = 0.0f;
    out.v4.y = 0.0f;
    out.v4.z = -(2.0f * far * near) / (far - near);
    out.v4.w = 0.0f;

    return out;
  }

  mat4 lookatf(vec3 eye, vec3 center, vec3 up)
  {
    static mat4 out;

    vec3 front, side, top;

    // relative to camera pose
    front = normalize(center - eye);
    side  = normalize(cross(front, up));
    top   = cross(side, front);

    out.v1.x = side.x;
    out.v1.y = top.x;
    out.v1.z = -front.x;
    out.v1.w = 0.0f;

    out.v2.x = side.y;
    out.v2.y = top.y;
    out.v2.z = -front.y;
    out.v2.w = 0.0f;

    out.v3.x = side.z;
    out.v3.y = top.z;
    out.v3.z = -front.z;
    out.v3.w = 0.0f;

    out.v4.x = 0.0f;
    out.v4.y = 0.0f;
    out.v4.z = 0.0f;
    out.v4.w = 1.0f;

    out = translate(out, -eye.x, -eye.y, -eye.z);

    return out;
  }

  dmat4 lookatd(dvec3 eye, dvec3 center, dvec3 up)
  {
    static dmat4 out;

    dvec3 front, side, top;

    // relative to camera pose
    front = normalize(center - eye);
    side  = normalize(cross(front, up));
    top   = cross(side, front);

    out.v1.x = side.x;
    out.v1.y = top.x;
    out.v1.z = -front.x;
    out.v1.w = 0.0f;

    out.v2.x = side.y;
    out.v2.y = top.y;
    out.v2.z = -front.y;
    out.v2.w = 0.0f;

    out.v3.x = side.z;
    out.v3.y = top.z;
    out.v3.z = -front.z;
    out.v3.w = 0.0f;

    out.v4.x = 0.0f;
    out.v4.y = 0.0f;
    out.v4.z = 0.0f;
    out.v4.w = 1.0f;

    out = translate(out, -eye.x, -eye.y, -eye.z);

    return out;
  }
}

